#include "def.h"
#include <stdio.h>
#include <cutil.h>
#include <cuda_runtime.h>
#include <cutil_inline.h>

typedef unsigned int uint;

/*
__global__ void binaryThreshold(uchar4 *p,uchar4 *o, int width, int height) {
  int x=blockIdx.x;
  int y=threadIdx.x;
  int z=y*width+x;
  int r=p[z].x, g=p[z].y, b=p[z].z;

  #define THRESHOLD 49152  
  if (r*r+g*g+b*b<THRESHOLD) { 
    r=g=b=255;  
  } else {
    r=g=b=0;
  }

  o[z].x=r;
  o[z].y=g;
  o[z].z=b;
  o[z].w=p[z].w;
}


__global__ void laplance(uchar4 *p,uchar4 *o, int width, int height) {
  int x=blockIdx.x;
  int y=threadIdx.x;
  int z=y*width+x;
  int r=p[z].x, g=p[z].y, b=p[z].z;
  float a[4];
  
  #define c(x,y,d) (p[(y)*width+x].d - p[z].d)
  #define A(x) (x*x)
 
  #define N(d) ( y>0 ? c(x,y-1,d) : 0 )
  #define W(d) ( x>0 ? c(x-1,y,d) : 0 )
  #define S(d) ( y<height-1 ? c(x,y+1,d) : 0 )
  #define E(d) ( x<width-1  ? c(x+1,y,d) : 0 )

  #define Cx(d) (W(d)>E(d)? W(d)-E(d) : E(d)-W(d))
  #define L(d) ((N(d) + W(d)) + S(d) + E(d))
  
  o[z].x=r + L(x)/4;
  o[z].y=g + L(y)/4;
  o[z].z=b + L(z)/4;
  o[z].w=p[z].w;
}

*/
__global__ void aniDiff(uchar4 *p,uchar4 *o,int width, int height, int K) {

  int x=blockIdx.x*blockDim.x + threadIdx.x;
  int y=blockIdx.y*blockDim.y + threadIdx.y;
//  int x=blockIdx.x, y=blockIdx.y;
  if (x>=width || y>=height) return;
  int i,z=y*width + x, m=threadIdx.z;

  float a[4],c[4];

  #define c(x,y,d) (p[(y)*width+x].d - p[z].d)
  #define sqr(x) (x*x)
 
  #define N(d) ( y>0 ? c(x,y-1,d) : 0 )
  #define W(d) ( x>0 ? c(x-1,y,d) : 0 )
  #define S(d) ( y<height-1 ? c(x,y+1,d) : 0 )
  #define E(d) ( x<width-1  ? c(x+1,y,d) : 0 )

  #define L(d) a[0]=N(d); a[1]=W(d); a[2]=S(d); a[3]=E(d); 
  
  if (m==0) { L(x) } else if (m==1) { L(y) } else { L(z) }
  
  #define Abs(x) (x>0 ? x : -x)
  for (i=0; i<4; i++) {
    c[i] = a[i];//-a[(i+2) % 4];
    c[i] = Abs(c[i])/K;
    //c[i]= 1/(1+A(c[i]));
    c[i]= exp(-sqr(c[i]));
  }

  float r=0;
  for (i=0; i<4; i++) 
    r+=a[i]*c[i];

  #define O(d) o[z].d=p[z].d + r/4;
//  #define O(d) if (p[z].d<150) o[z].d=0; else o[z].d=255;
  if (m==0) { O(x) } else if (m==1) { O(y) } else { O(z) }
  o[z].w=p[z].w;
}

// implement following functions
//
bool kernel_2d_convolution_gpu(uchar4 *_in, uchar4 *_out, int width, int height,int K, int SIZE)
{
  //cudaGetDeviceProperties();
  int N=width*height;
  uchar4 *p,*o;
  int size=sizeof(uchar4)*N;

  cudaMalloc((void**)&p,size); 
  cudaMalloc((void**)&o,size);
  
  cudaMemcpy( p, _in, size, cudaMemcpyHostToDevice);
  
  //binaryThreshold<<<width,height>>> ( p, o, width, height);
  //laplance<<<width,height>>> ( p, o, width, height);

  dim3 gridSize, blockSize;
//  gridSize.x=width; gridSize.y=height; gridSize.z=1;
//  blockSize.x=3; blockSize.y=blockSize.z=1;
  gridSize.x=(width-1)/SIZE+1;  gridSize.y=(height-1)/SIZE+1; gridSize.z=1;
  blockSize.x=SIZE;             blockSize.y=SIZE;             blockSize.z=3;
//  printf("\nGrid: %d %d %d ",gridSize.x,gridSize.y,gridSize.z);
//  printf("Bloc: %d %d %d\n",blockSize.x,blockSize.y,blockSize.z);
  aniDiff<<<gridSize,blockSize>>> ( p, o,width,height, K);

  cudaMemcpy( _out, o, size, cudaMemcpyDeviceToHost);

  cudaFree(p);
  cudaFree(o);
	return true;
}
